Homework 3
Int. Data Analytics and Machine Learning (CMDA/CS/STAT 4654)
Catherine Squillante
Instructions
This homework covers the multiple linear model and stepwise model selection lectures. All work must be submitted electronically. For full credit you must show all of your steps. Use of computational tools (e.g., R) is encouraged; and when you do, code inputs and outputs must be shown in-line (not as an appendix) and be accompanied by plain English that briefly explains what the code is doing.
Problem 1: Infant nutrition (20 pts)
This question involves data from a study on the nutrition of infants and preschool children in the north central region of the United States of America It is available on the course web page as nutrition.csv, and contains 72 observations of boys’ weight/height ratio (woh) for equally spaced values of age in months.
- (3 pts) Plot the data, and overlay the least squares line and a 95% prediction interval in the range of the data. Comment on the goodness of fit.`
nutrition <- read.csv("~/Documents/Documents/Machine Learning/nutrition.csv")
fit <- lm(woh ~ age, data=nutrition)
p <- predict(fit, newdata = nutrition, interval = "prediction")
plot(nutrition$age, nutrition$woh, main="prediction interval", xlab = "age", ylab="woh")
x <- seq(0,max(nutrition$age),length=length(nutrition$age))
lines(x, p[,1])
lines(x, p[,2], col=2, lty=2)
lines(x, p[,3], col=2, lty=2)
The fit tells us that there is a 95% chance new data points will fall in between the prediction interval (represented by the red lines on the graph.) The r squared value is .8199, which means approximately 82% of the variability of the response data around the mean is explained by the model.
- (7 pts) The authors of the study have reason to believe that the observations fall into two groups: (1) the first seven boys and (2) the remaining 65. By introducing an appropriate dummy variable and interaction term, find the least squares equations of these lines, and give the 95% confidence intervals of the corresponding intercepts and slopes. Are these significantly different from zero?
cat <-nutrition$age < 6.5
early <- as.factor(cat)
reg <- lm(nutrition$woh ~ nutrition$age + early + nutrition$age:early)
conf <- confint(reg, level = .95)
##intercept, slope, CI for first group:
c(reg$coefficients[1] + reg$coefficients[3], reg$coefficients[2] + reg$coefficients[4])## (Intercept) nutrition$age
## 0.43304762 0.04342857##int CI
c(conf[1]+conf[3], conf[1,2]+conf[3,2])## [1] 0.3756907 0.4904045##slope CI
c(conf[2]+conf[4], conf[2,2]+conf[4,2])## [1] 0.03120052 0.05565662##intercept, slope, CI for second group:
c(reg$coefficients[1], reg$coefficients[2])## (Intercept) nutrition$age
## 0.719655995 0.004157843##int CI
c(conf[1], conf[1,2])## [1] 0.7056933 0.7336187##slope CI
c(conf[2], conf[2,2])## [1] 0.003836169 0.004479518Yes, these are statistically different from 0 because all of the values are statistically significant based on the lm results.
- (5 pts) Overlay the two new lines and corresponding prediction intervals on your plot from part (a) in such a way as they cover only their respective
ageranges (i.e., so that they do not overlap). Comment on the goodness of fit of this new model.
plot(nutrition$age, nutrition$woh, main="prediction interval", xlab = "age", ylab="woh")
x <- seq(6.5,max(nutrition$age),length=length(nutrition$age))
i1 = reg$coefficients[1] + reg$coefficients[3]
s1 = reg$coefficients[2] + reg$coefficients[4]
i2 = reg$coefficients[1]
s2 = reg$coefficients[2]
ci1 = conf[1]+conf[3]
si1 = conf[2]+conf[4]
ci2 = conf[1,2]+conf[3,2]
si2 = conf[2,2]+conf[4,2]
segments(x0 = -1, y0=i1 + -1*s1, x1=6.5, y1=i1+6.5*s1)
segments(x0 = 6.5, y0=i2 + 6.5*s2, x1=80, y1=i2+80*s2)
segments(x0 = 6.5, y0=conf[1] + 6.5*conf[2], x1=80, y1=conf[1]+80*conf[2], col=2, lty=2)
segments(x0 = 6.5, y0=conf[1,2] + 6.5*conf[2,2], x1=80, y1=conf[1,2]+80*conf[2,2], col=2, lty=2)
segments(x0 = -1, y0=ci1 + -1*si1, x1=6.5, y1=ci1+6.5*si1, col=2, lty=2)
segments(x0 = -1, y0=ci2 + -1*si2, x1=6.5, y1=ci2+6.5*si2, col=2, lty=2)
- (5 pts) Explain a test which allows us to choose between the fits in (a) and (b). Perform the test and explain the result.
fit = lm(nutrition$woh ~ nutrition$age)
anova(fit, reg)## Analysis of Variance Table
##
## Model 1: nutrition$woh ~ nutrition$age
## Model 2: nutrition$woh ~ nutrition$age + early + nutrition$age:early
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 70 0.188623
## 2 68 0.042335 2 0.14629 117.49 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform the partial F test. This tells us whether the added predictors are worthwile. The F value in obtained from the anova test above is very statistically significant. The null hypothesis is rejected, meaning the added predictor is useful and the fit from part(b) is better.
Problem 2: Beef – its what’s for dinner (15 pts)
In 1988, US cattle producers voted on whether or not to each pay a dollar per head towards the marketing campaigns of the American Beef Council. At the time of this vote, the council’s TV campaign featured a voice-over by actor Robert Mitchum, using the theme “Beef – it’s what’s for dinner.” To understand the vote results (it passed), the Montana state cattlemen’s association looked at the effect of the physical size of the farm and the value of the farms’ gross revenue on voter preference. The data (in beef.csv) consist of the vote results (%YES), average size of farm (hundreds of acres), and average value of products sold annually by each farm (in K$) for each of Montana’s 56 counties.
- (5 pts) Plot the data and comment on what you see. How will this effect your analysis?
beef <- read.csv("~/Downloads/beef.csv")
plot(beef)
This plot shows that there is a pretty strong positive correlation between size and val.There is a moderate negative correlation between yes and size. And a weak negative correlation between yes and val. There is multicollinearity because the predictors (size and val) are correlated with themselves. They are actually more correlated with themselves, so it will be harder to make confident inference about these predictors. There is also some bunching (more data points for val at the lower x values), so we will probably want to take the logarithm.
- (5 pts) Fit a regression model for
YESwith bothSIZEand \(\log(\)VAL\()\) as covariates. Interpret the results. What regression assumptions might we have violated here?
attach(beef)
reg = lm(YES ~ SIZE + log(VAL))
summary(reg)##
## Call:
## lm(formula = YES ~ SIZE + log(VAL))
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.2200 -5.3463 0.5873 7.7057 28.0153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.71445 9.95022 9.418 6.55e-13 ***
## SIZE -0.24275 0.08743 -2.777 0.00758 **
## log(VAL) -2.38614 2.79193 -0.855 0.39659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.36 on 53 degrees of freedom
## Multiple R-squared: 0.2606, Adjusted R-squared: 0.2327
## F-statistic: 9.342 on 2 and 53 DF, p-value: 0.0003346SIZE is a significant predictor, as summary shows the low p value. However, VAL may not be a good predictor since is has a high p value and is not very significant. This violates the assumption that the predictor variables are independent from one another leading to large standard errors.
- (5 pts) Find a better model: does the effect of
SIZEchange depending on \(\log(\)VAL\()\)? What is your estimate of the effect onYESof a unit changeSIZE? Interpret your conclusion.
attach(beef)## The following objects are masked from beef (pos = 3):
##
## SIZE, VAL, YESreg1 = lm(YES ~ SIZE*log(VAL))
summary(reg)##
## Call:
## lm(formula = YES ~ SIZE + log(VAL))
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.2200 -5.3463 0.5873 7.7057 28.0153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.71445 9.95022 9.418 6.55e-13 ***
## SIZE -0.24275 0.08743 -2.777 0.00758 **
## log(VAL) -2.38614 2.79193 -0.855 0.39659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.36 on 53 degrees of freedom
## Multiple R-squared: 0.2606, Adjusted R-squared: 0.2327
## F-statistic: 9.342 on 2 and 53 DF, p-value: 0.0003346summary(lm(YES ~ SIZE))##
## Call:
## lm(formula = YES ~ SIZE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.8015 -5.8995 0.5392 7.9553 28.4878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 85.48991 2.52366 33.875 < 2e-16 ***
## SIZE -0.28940 0.06813 -4.248 8.57e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.33 on 54 degrees of freedom
## Multiple R-squared: 0.2505, Adjusted R-squared: 0.2366
## F-statistic: 18.04 on 1 and 54 DF, p-value: 8.572e-05Yes, the effect of size changes depending on log(VAL). The summary shows that the SIZE: log(VAL) relationship is significant based on its small p value. Based on the second R output above, there is a negative correlation between size and yes. For every unit increase in size, there is a ~.289 decrease in yes.
Problem 3: Weighted regression (20 pts)
Consider a situation where each \((x_i, y_i)\) pair is given a weight \(w_i > 0\), for \(i=1,\dots,n\) measuring its importance relative to the other pairs in the data. Otherwise the multiple linear regression setup is the same as always: \(x_i\) and \(\beta\) are \(p\)-vectors, the variance is constant, etc.
(5 pts) Ignoring for the moment the typical Gaussianity assumption on the noise, consider the following loss function similar as a weighted analog of the ordinary least squares. \[ \sum_{i=1}^n w_i (y_i - x_i^\top \beta)^2 \] Find an expression for the setting \(\hat{\beta}\) that minimizes that loss. Hint: You might find it helpful to vectorize the loss function above.
(4 pts) Now, suppose that your data was comprised of \(n\) unique \(p\)-variate \(x_i\) predictors, and at each one of those settings you had \(w_i \geq 1\) observations \(y_{i1}, \dots, y_{iw_i}\). These are \(w_i\) replicates observed at \(x_i\). Write down the full multiple linear regression model for such data, assuming \(\epsilon_{ij} \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0, \sigma^2)\), and derive an expression for the corresponding the log likelihood. Hint: consider using a double sum.
(5 pts) Show that the MLE \(\hat{\beta}\) from the likelihood in (b) may be found by applying your results from (a–b) using \(y_i \equiv \bar{y}_i = \frac{1}{w_i} \sum_{j=1}^{w_i} y_{ij}\). That is, although you have \(\sum_{i=1}^n w_i \equiv N\) data points, and potentially \(N \gg n\), once you have the \(\bar{y}_i\)’s your calculations need only involve working with \(n\) quantities, and you’ll get the same answer as if you had performed the calculations with all \(N\) \(y\)-values.
(6 pts) However, show that using the \(n\) averages \(\bar{y}_i\) is not sufficient for estimating \(\hat{\sigma}^2\). That is, you won’t get the same answer using just the \(\bar{y}_i\)’s as you would by working with all \(N\) quantities. What must you additionally know about each collection of replicates \(y_{i1}, \dots, y_{i w_i}\) so that your calculations involve working with \(2n\) quantities rather than \(N\)? Derive \(\hat{\sigma}^2\) using those \(2n\) quantities.
Problem 4: Cook’s distance (20 pts)
Consider the linear model \(Y = X\beta + \epsilon\), where \(X\) is a known \(n \times p\) matrix, \(\epsilon \sim N_n(0,\sigma^2I)\), and let \(x_i^\top\) denote the \(i\)th row of \(X\).
- (5 pts) Let \(X_{(-i)}\) denote the \((n-1) \times p\) matrix obtained by deleting the \(i^\mathrm{th}\) row of \(X\), and suppose that this matrix is also of full rank \(p\). Show that \[
X_{(-i)}^\top X_{(-i)} = X^\top X - x_ix_i^\top.
\]
(7 pts) Let \(\hat{\beta} = (X^\top X)^{-1}X^\top Y\) denote the MLE in the model with all \(n\) observations, and let \(\hat{\beta}_{(-i)} = (X_{(-i)}^\top X_{(-i)})^{-1}X_{(-i)}^\top Y_{(-i)}\) denote the MLE in the model with the \(i^\mathrm{th}\) observation deleted. Show that \[ \hat{\beta} - \hat{\beta}_{(-i)} = \frac{1}{1-h_i}(X^\top X)^{-1}x_i e_i, \] where \(h_i\) is the leverage of the \(i^\mathrm{th}\) data point and \(e_i\) is our usual residual. Hint: you might find the Sherman–Morrison–Woodbury formula useful.
(6 pts) Cook’s distance of the observation \((x_i,y_i)\) is defined as \[ D_i = \frac{1}{ps^2}(\hat{\beta}_{(-i)} - \hat{\beta})^\top X^\top X(\hat{\beta}_{(-i)} - \hat{\beta}), \quad i=1,\ldots,n, \] where \(s^2 = \|Y - X\hat{\beta}\|^2/(n-p)\). It may appear that we need to fit \(n+1\) linear models in order to calculate all of the Cook’s distances. Deduce, however, that \[ D_i = \frac{1}{p}\Bigl(\frac{h_i}{1-h_i}\Bigr)\tilde{r}_i^2, \] where \(\tilde{r}_i = \frac{e_i}{s\sqrt{1-h_i}}\) is (very similar to) the \(i^\mathrm{th}\) studentized fitted residual (given in class). Actually, this is called the internally studentized residual, whereas our version in class is the externally studentized one.
